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Abstract 

Mass and spin are the only two parameters needed to completely characterize black 
holes in General Relativity. However, the interaction between black holes and their 
environment is where complexity lies, as the relevant physical processes occur over a 
large range of scales. That is particularly relevant in the case of super-massive black 
holes (SMBHs), hosted in galaxy centers, and surrounded by swirling gas and var- 
ious generations of stars. These compete with the SMBH for gas consumption and 
affect both dynamics and thermodynamics of the gas itself. How the behavior of such 
fiery environment influence the angular momentum of the gas accreted onto SMBHs, 
and, hence, black-hole spins is uncertain. We explore the interaction between SMBHs 
and their environment via first 3D sub-parsec resolution simulations (ranging from 



~ 0.1 pc to ~ 1 kpc scales) that study the evolution of the SMBH spin by including 
the effects of star formation, stellar feedback, radiative transfer, and metal pollution 
according to the proper stellar yields and lifetimes. This approach is crucial to in- 
vestigate the impact of star formation processes and feedback effects on the angular 
momentum of the material that could accrete on the central hole. We find that star for- 
mation and feedback mechanisms can locally inject significant amounts of entropy in 
the surrounding medium, and impact on the inflow inclination angles and Eddington 
fractions. As a consequence, the resulting trends show upper-intermediate equilibrium 
values for the spin parameter of a ~ 0.6 — 0.9, corresponding to radiative efficiencies 
e :i^ 9% — 15%. These results suggest that star formation feedback taking place in the 
circum-nuclear disk during the infall cannot induce alone very strong chaotic trends in 
the gas flow, quite independently from the different numerical parameters. 
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Introduction 



Black holes (BH) in Einstein's theory of general relativity are described by only three parameters: 
mass, spin and charge. Astrophysical black holes are even simpler, as charge can be neglected. So, 
besides their masses, Mbh, they are completely characterized by their dimensionless spin parame- 
ter, a = c Jbh/G Mgjj, where c is the speed of light, G the gravi t ation constan t , Jrh is the angular 
momentum of the black hole, and < a 5, 1 (e-g - iKern 1 19631 : IPetersI 1 19641 : iBoyer & Lindquist 



19671 : ICartedll968l: lFelicelll968l: iBardeen et al.lll972l) . 

Many previous theoretical effo rts have focused on the mass growth of SMBHs that reside in galaxy 



centers and power quasars (e.g. Lvnden-Belllll969l: 



19961: lOhosh & Abramowicz 



LivioetaU1999: 



1997 



Pi Matteo et al.ll2005 



Blandford & Znaieklll977l:lAbramowicz et al.lll978 



Bardeen 



Moderski et al 



2008 



Bower et al 



2006l.l2007l:lsiiacki et al.ll2007l:lBooth & Schavell2009l:lsiiacki et al 



1970l:lshakura & SunvaevI 



Rees et 



1998; 



all 19821: 



Wilson & Colbert 



Silk & ReesI 



2006 



19981: 



973, 



1995 



Haehnelt et al, 



Croton et al. 



2006 



Li et al 



1976; 



Moderski & Sikora 



19981 : 



20061: 



200511201 ll : lHopkins & OuataertlboiO 



Fabjan et al 



20 lOt Dubois et al .120121 and references therein) 



20091 : iHopkins et al. 



Thus, we will focus here on the processes taking place in their environments and on the consequent 

implications for the other property of SMBHs: spin, whi ch is of no less importanc e. 

Indee d, the spin of a black hole affects the efficiency (e.g. lNovikov & Thomdll973l:lBambi & Barausse 
20 111 and references therein) of 'classical' accretion processes. In radiatively efficient thi n ac- 
cretion disks the mass-to-energy conversion efficiency, e, equals (|Novikov & Thomelll973b e = 
1 — E/c^, where E is the binding energy per unit mass of a particle at the last stable orbit. The 
closer the last stable orbit is to the horizon, the higher the mass-to-energy conversion efficiency. 



which increases from e ^ 5.7%, for a non-rotating (a = 0) hole, to e ~ 42%, for its maximally 
rotating counterpart (see figure [1}. 



The mass-to-energy conversion directly af fects the mass growth rate of black holes (|Soltanlll982 : 
Yu & Tremainell2002l: IShankar et al.ll2009|) according to 



M 
t 



/] 



>Edd 



Edd 



(1) 



where /sdd is the Eddington efficiency, and L^dd the Eddington luminosity: high radiative effi- 
ciency, e, implies slow growth, as more mass is radiated away instead of being fed into the black 
hole. For instance, for a black hole accreting at a fraction f^dd of the Eddington rate, ^Edd, assum- 
ing constant e, the mass increases with time as: 



M(t) 



M(0)exp(-^^^'^^^"') 



t 






tEdd 



(2) 



where M(0) is the initial mass, and the Eddington timescale, tsdd, is c^ divided by the Eddington 
luminosity per unit mass, 

tEdd = ^^—^ 0.45 Gyr, (3) 



AnGrrip 

with ar Thomson cross section, and rrip proton mass. Moreover, the higher the spin parameter, 
a, the higher e, implying longer timescales to grow the black-hole mass (see also later figure |2]). 
Going from e ~ 5.7% to e ^ 42%, the difference in mass amounts to 6 orders of magnitude at 
t = t^dd, for fEdd = 1- The typical spin therefore impacts the overall mass growth of SMBHs. 
This is particularly imp ortant for the case of SMBHs powering the highest-redshift (z > 6) quasars 
king & I 
lion year. 

Although the initial mass of the SMBH plays a crucial role in this game, all the SMBH formation 
models proposed to date require significant and fast mass accretion to result in lO^Mo SMBHs 



mr)( 

m. 



(iKing & Pringldl2007l) . since these SMBHs must grow in mass by gas accretion in less than a bil- 



Fig. 1. — Radiative efficiency as a function of the spin parameter, a (IBardeen et al.l 1 19721 : 



Novikov & Thornelll973|) . Positive values of a correspond to SMBHs co-rotating with their accretion disks, 
while negative values for a indicate counter-rotating SMBHs. 



at redshift > 6. As an example, under the assumption of a constant spin and radiative efficient 
accretion, a non-rotating (a = 0, e ~ 5.7%) SMBH accreting at its Eddington limit (/Edd = 1) can 
grow its mass of 6 orders of magnitude in only < 0.4 Gyr. A maximally rotating (a = 1, e ~ 42%) 
SMBH would accrete (at /sdd = 1) the same mass in more than a billion year. Note however that 
assuming a constant a = is highly unphysical, since non rotating SMBHs are spun-up by every 
accretion event with finite angular momentum. Fortunately, the constraint on the spin parameter 
are not too s tringent, since (see figure]!]) the radiative effi ciency is a very steep function of a only 



for a > 0.9 (|Novikov & Thomelll973l: iKing et al.ll2008al) . with typical values varying from a few 
per cents up to e si 10% for a < 0.7, and slightly increasing to e ~ 15% for a ^ 0.9. As an 
example, from equation ^ one can estimate that a SMBH with a Ri 0.8 (e ^ 11.5%) would grow 
by 6 orders of magnitude in roughly two times the time needed by a non-rotating SMBH still well 
within the age of the Universe at 2; ?» 6 (of about 1 Gyr). 

The magnitude and orientation of SMBH spins also affect the frequency of SMBHs in galaxies, 
via 'gravitational recoil' . When two black holes coalesce, they may recoil due to anisotropic emis- 
sion of gravitational waves. If this recoil were sufficiently violent, the newly merged hole would 
escape from the host galaxy. Recent breakthroughs in numerical relativity have allowed reliable 
computations of recoil velocities. Non-spinning black holes are expected to recoil with velocities 
below 200 km s~^. The recoil is m uch larger, up to thousands km s~^, for SMB Hs with large spins 



in non-aligned configurations (e.g. ICampanelli et al.ll2007l: iLousto et al 



Final ly, the value of the spin parameter, a, in a Kerr black hole (|Kenill963l : iBoyer & Lindquist 



2012) 



967h also determines how much energy is in princ iple extractable from the black hole itself 



2005h in relation to jet production (but see 



(iBlandford & Znaieklll977l:lBlandford & Pavndl 19821). a nd it may be important (as argued by e.g. 
MacDonald & Thomelll982l: iModerski & Sikoralll996[ lohosh & Abramowic3ll997l: iMcKinnev 



Livio et al. 



19991 



i 



In general, since SMBH growth i s supposed to be led by mass accretion (|Soltanlll982l : lYu & Tremaine 
2OO2I : Ishankar et al.ll2009l boilb . it is the way SMBHs accrete gas that will determine their spin. 
Spin-up is a natural consequence of prolonged disk-mode accretion: any hole that has, for instance, 
doubled its mass by capturing material with constant angular momen tum axis (cohe rent accretion) 
would end up spinning rapidly, close to the maximum allowed value (|Thomelll974|) . For example, 
if a BH accretes at values of order ~ 10% the Eddington limit (i.e. /sdd = 0.1), then it can be 
spun up to high values in less than 10 Myr (figure |2]). The radiative efficiency evolves much more 
slowly, since as discussed previously, it is a shallow function of the spin until a ~ 0.9. 
On the other hand, a black hole which had gained its mass from capturing many low-mass objects 
in randomly-oriented orbits or from small (an d short) uncorrelated accreti on episodes (chaotic ac- 
cretion), would keep small spin (a < 0.2, e.g. lKing & Pringlell2006ll2007l) . 

The main unknown here is therefore the typical angular momentum of the material feeding the 
SMBH: does the vector maintain a roughly constant preferential direction, or does the direction 
change chaotically? 



Since SMBHs are expected to have grown in mass mostly during active phases as quasars, the evo- 
lution of SMBH spins has to be addressed by simulating the typical environmental conditions ex- 
pected around quasars, where the properties of the acc retion flow are established. When accretion 
is triggered by galaxy mergers, as expected for quasars (jSanders & Mirabell 19961 lOownes & Solomon 
19981 : bavies et allboOBJ I2OO4I : bownes & EckartlboOTi) . the material that feeds a SMBH is ex- 
pected to assemble into a circum-nuclear disk (CND), as also found by numerical simulations (e.g. 
Mayer et al.ll2007h . These are therefore the most relevant circumstances to be explored. CNDs 
are also associated to very high levels of star formation that are expected to influence SMBH spin 
evolution: local star formation injects energy and momentum into the gas near a SMBH, possibly 
breaking the coherency of the gas inflow. A quasar could be fueled by a sequence of well-separated 
short -lived events with randomly distributed angular momentum vectors, leading to a chaotic gas 
flow king et alJl2008bl) . 

In this work we will try to understand if a coherent circum-nuclear disk around a SMBH can be 
affected by the gaseous and stellar processes taking place during infall. Moreover, we will try 
to address if these phenomena can break or alter the coherency of the matter flow, independently 
from larger- scale ef fects, like mergers or SMBH-SMBH interactions (as already studied by e.g. 
Hopkins et aDbo lib . 

Therefore, we simulate the gaseous and stellar environment around a 4 x 10^ Mq SMBH placed 
in a circum-nuclear disk of 10^ Mq, by using the parallel tree/SPH code Gadget-3, a n improved 
version of the publicly available code Gadget- 2 dSpringel et al.ll200ll: ISpringellbOlOl) . including 



cooling (IMaio et al 



2OO7L l2010aL I2OI ih. star formation dSpringel & Hemquistll2003ah . stellar, and 



radiative feedback (|Petkova & SpringellbOllh . We make sure that our spatial resolution (from 0.1 



to 1 pc) is sufficient to resolve the 'sphere of influence' within which the SMBH dominates gravity 
(^ 10-20 pc), allo\v ing us to study the dynamics of the inflowing material (e.g. iDottiet al.ll2010 : 
Hopkins et al.ll201 ih . This also allows us to model gas cooling, star formation, and stellar feedback 
accurately. In our investigation we focus on the effect of local perturbatio ns, neglecting large- scale 
instabilities that can alter the angular momentum of the nuclear gas (e.g. iHopkins et al.ll201 ih . 

The paper is organized as follows: in section |2l we introduce and describe the simulations we 
performed; in section l3l and 13.31 we outline the technique used to address spin evolution and show 
our main results; in section IH we comment and conclude. 



2. Simulations 



In the following we will describe the code used, the physical processes included in our treatment, 
and the set of 3D simulations performed. 



2.1. Implementation 



We run 3D N-body/hydro simulations by usin g the parallel tree, SPH code Gadge t3, an extended 
version of the publicly available code Gadget2 (|Springel et alJl2005uSpringelll2010). Besides grav - 



ity and hydrodynamics, the code includes gas cool ing down to low temperatures (IMaio et al. 



2003h. inspired by the works of iKatz & GunnI (1199 ih : 



2010ah and a multi-phase model f or star formation (ISpringel & Hemguist 2003bl:lHemauist & Springel 



Cen & Ostrikei 



Katz et alj (1 19921 . 1 1 9961) . Ad ditionally, mechanical feedbac k (|Aguirre et al 



(SN) is taken into account (ISpringel & Hemguist 



(RT) prescriptions (IPetkova & Springelll2009l 12011 



2007 



1992c 



aib|): iKatzl 



iim); 



I2OOII) from supemovae 



2003bh. We also implern ent radiative transfer 



Petkova & Maid 12012 ) to model the effects 



onto the surrounding medium of photo-ionization and photo-heating from stellar radiation and 
gaseous radiative processes. 

The main agents of gas cooling are atomic resonant transitions that can rapidly bring the temper- 
atures down to ~ 10"* K. In me tal-rich enyironrnents, at omic fine- structure transitions can lower 
temperatures further to ~ 10 K dMaio et al.ll2007l . l2010bh . We will additionally show that detailed 
metallicity evolution can slightly influence the resulting star formation and feedback processes, 
but has no drastic effects on the spin evolution of the central SMBH and helps stabilizing the disk 
against star- formation-induced chaotic motions. 

When the medium cools and gets denser than a given density threshold of 10^ cm~^, stochastic 
star formation is assumed, cold gas is gradually converted into stars, and entropy feedback is in- 
jected into the surrounding environment, leading to a self-regulated star forming regime. The only 
free parameter in this model is the maximum star formation timescale, t°, which determines the 
amplitude of the typical star formation timescale, U, at any gas density, p, above the mass-density 
threshold, pth. 



t^P) = tU- 



p_ 
Pth 



-1/2 



(4) 



In some astrophysical problems, it is possible to fix t° by matching with observations!!!. In strongly 
dynamically evolving systems, such as intense star-bursts in galactic nuclear regions, however, it is 
impossible to define a unique characteristic timescale for star formation. Thus, we will assume two 
extreme ti mescales, a very long o ne (t° = 300 Myr), and a very short one (t° = 1 Myr), as limiting 
cases (e.g. !Krumholz et al.l!201 ih . This means that the maximum timescales at the threshold for 
star formation are 300 Myr and 1 Myr, respectively, while at p > pth, the time over which stars can 
form, t*(p), becomes shorter and shorter. Moreover, we will also explore the implications of inter- 
mediate, more realistic values (t° = 6, 10, 20, 50, 100 Myr) in terms of star formation efficiency, 
and the effects of full stellar evolution with consistent metal enrichment calculations. 



' For example, in galact ic disks, i^ = 2.1 Gyr allows us to fit the Kennicutt-Schmidt law (lKennicutllll998l : 
Springel & HernquisjEooib!) . 



At each timestep, At ^ 1 Myr, coUisionles star particles are spawned with a probability distribu- 
tion 

P. = — \l- e-(^-^)^^*/**l , (5) 



where m is the current mass of the gas particle, m^ is the mass of each star particle, /3 ~ 0.1 is 
the fractioncl of short-lived stars which die in stantly as supernovae, and x is the cold-cloud mass 
fraction, ranging roughly from 0.85 up to ~ 1 (|Kennicuttlll998l : ISpringel & Hemquistil2003bl) . The 
star mass is given by m^ = rrio/Ng, with A^^ an integer number in dicating how many stars each 
gas particle can generate - usually of the order of unity or a few (jSpringel & Hemquistll2003b : 
Hemquist & Springelll2003h ) - and rriQ the initial gas particle mass, which is decreased (and re- 
duced to the current gas mass m) accordingly to the spawned stars. The resulting star formation 
rate is simply given by M* = (1 — /3)xm/ti,. The initial position and velocity of each star particle 
are assumed to be the same as the parent gas particle. We note that relations dH) and © imply 
higher probability of forming stars at higher densities. Once the star forming regime kicks in, su- 
pernova energy heats the ambient medium and cold clouds are evaporated by thermal conduction: 
this process is parameterized via the evaporation factor, Aq, dependent on the environment. As- 
suming that thermal instabilities in the gas become active at ~ 10^ K and typical SN temperatures 
are Tsn ~ 10^ K, the evaporation factor is Aq = Tsat/IO^K ~ 10^. However, different SN and 
cooling feature s would give different values, t hus, we will assume Aq = 10^, lO'', and 10^. For fur- 
ther details see ISpringel & HemquistI (|2003b|) . Similarly, outflows are phenomenologically taken 
into account according to the following probability distribution: 



Pv 



1 -e" 



-7?{l-/3)a;Ai/i* 



(6) 



where the wind efficiency, 77 ~ 2 (|Martinll 19991 : ISpringel & Hemquistll2003bh . is the ratio between 
the mass-loss rate and the star formation rate of each gas particle. Wi nds carry an additional ki- 
netic energy that is a fraction, x — 0.25 (|Springel & HemquistI l2003b|) . of the supernova energy. 
Therefore, the wind velocity, v^j, is v^j = \/2xe^T^/rj , where ^sn = 4 x lO'^^ erg Mq^ is the 
supe rnova energy injected in the rn edium for each solar mass in stars formed (assuming a Salpeter 
IMF, ISpringel & HemquistI l2003b|) . The corresponding wind velocity has typical magnitudes of a 
few hundreds km/s in a direction that can be randomly chosen either on the unit sphere (isotropic 
winds) or along the axis uniquely defined by v x V 0, with velocity vector v and gravitational 
potential </> (axial winds). 

The photo-ionization and photo-heating of the gas is followed by a radiative transfer scheme ex- 
tended to a multi -frequency regime, with four frequency bins in the range 13.6 - 60 eV, and cor- 
responding to the ionization potential of neutral hydrogen (13.6 eV), neutral helium (24.6 eV), 



^Here a Salpeter IMF has been assumed over the range [0.1,100] M^ 



singly ionized helium (54.4 eV), and one arbitrarily higher value (60 eV). The radiative trans- 
fer equations are coupled to hydrodynamics and solved self-consistent ly taking account both gas 



cooling and radiative heating from stellar sources (for further details see lPetkova & Springelll2009 



2011 



Petkova & Maid 1201 2|) . Ionizing radiation from stars and gas is diffused throughout the 



simulation volume, interacting with both hydrogen and helium. The stellar emissivity is assumed 
to have a black-body spectrum with effective temperature T^s = 5 x 10^ K (usual emission from 
short-lived massive stars). The gas emissivity due to cooling processes is also a black-body spec- 
trum with effe ctive temperature equal to its own temperature. 

A few studies (|Wada et al.ll2009l : ISchartmann et al.i2009U201 1|) attempted analyses of the gas prop- 
erties around SMBHs. None of them, though, was suited to study how star formation influences 
the evolution of SMBH spins, which is the main goal of our investigation. Previous works, in fact, 
made several simplifications on the physical processes considered, for instance neglecting either 
radiative transfer effects or the fundamental link between stellar evolution and the gaseous envi- 
ronment where stars form. Our simulations, instead, follow star formation and stellar evolution 
explicitly, so that SN explosions are driven from the hydrodynamical evolution of the gas. The 
ensuing feedback effects are therefore fully self-consistent with the dynamics and thermodynam- 
ics of the circum-nuclear disk. These features allow us to investigate accurately how stellar and 
SN feedback affect the angular momentum of the gas near SMBHs, and, as a consequence, how 
star-bursts impact SMBH spins. 



2.2. Initial conditions 



We initialize every simulation with a SMBH of 4 x lO^M© at rest, at the center of a purely gaseous 
CND. The SMBH is treated as an active evolving coUisionles particle, and reacts to the gravita- 
tional torques exerted by the stellar and gaseous background. 



The SMBH and the disk are embedded in a spherical stellar bulge (IMagorrian et al.lll998|) with 
total mass of 2 x lO^M© and scale length b = 10^ pc, always modelled using 10^ particles of 



2_2<10^Mq. The bulge follows a Hemquist mass density profile (lHemquistill990l : lTremaine et al 



1994|) . which reproduces the i?^/^-law well (|de Vaucouleurslll948h and mimics the potential of the 
host galaxy. 

The disk is assumed to have a mass of IO^Mq (25 times larger than the central SMBH) and is 
modelled using 2 x 10^, 2 x 10^, 2 x 10^, or 2 x 10^ gas particles in our low, standard, high, and 
very-high resolution runs, such that the corresponding gas particle masses are 5000 M©, 500 M©, 
50 Mq, and 5 Mq. The CND is initially composed only of gas particles, since we aim at modelling 
the evolution of SMBHs in a recently formed CND. We note, however, that the CND becomes 
populated of young stars as the simulation evolves. The absence of an old population of stars in 
our initial conditions does not affect our results: only young stars, through their powerful winds 



or their deaths as supemovae, can significantly alter the properties of the gas. The disk initially 
follows an exponential surface density profile, and has a scale length of 100 pc. It is pressure 
supported in the vertical direction, with an aspect ratio H/ R = 0.1. The disk internal energy is 
initially fixed according to a polytropic equation of state with adiabatic index 7 = 5/3 and cosmic 
composition for H and He. 

We evolve adiabatically our initial conditions for fti 100 Myr, to ensure that the system reaches 
a stable configuration without fragmenting nor developing any strong asymmetry and instability, 
before turning on cooling, star formation, winds, supernova and radiative feedback. 



2.3. Runs 



The runs we performed involve the study of several parameters and physical processes and we refer 
to Tabled] for the general outline of the features in the different simulations^. 
The first, second, and third blocks refer to the simulations with axial wind feedback and evap- 
oration factor Ao = 10^, Aq = 10'', and Aq = 10^, respectively. For each case we run four 
simulations: two with t° = 1 Myr (one with metallicity Z = 0, and H and He cooling down to 
~ 10^ K, and one with additional metal low-T cooling down to ~ 10 K for an assumed constant 
Z = Zq) and two with t° = 300 Myr (one with Z = 0, H and He cooling down to ~ lO'^ K, and 
one with additional metal low-T cooling down to --^ 10 K with Z = Zq). 

The corresponding complete, more detailed runs including (i) self-consistent metal production 
from stellar evolution according to the proper yields and lifetimes, (ii) exact abundance calcula- 
tions, and (iii) metal-dependent cooling in the range ~ 10 — 10^ K, will be also discussed in the 
following sections for the two limiting cases t° = 1 Myr and t° = 300 Myr (both with Aq = 10^). 
The resulting consequences from the cooling down to low temperatures on disk stability and frag- 
mentation are readily seen in the density maps of figure |3] 

In order to check the effects of outflows, we also run additional simulations with isotropic winds 
(IW), for the case with Aq = 10^, cooling down to ~ 10^ K, and either t° = 1 Myr or t° = 
300 Myr. 

We carry out convergence studies over the dependence on mass and spatial resolution by perform- 
ing four additional runs, with Aq = 10^, t° = 1 Myr and no low-T cooling. The "standard" runs 
have a gravitational softening length of 1 pc and 2 x 10^ gas particles, with corresponding gas 



^ In the runs, we used 32 SPH neighbours for the density calculations, however, we checked the bahaviours up to 
64 neighbours, as well, and found no significant difference in the results (mass profiles, star formation rates, angular 
momentum, relative angles, accretion rates, etc.). 
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Fig. 2. — Time evolution of the properties of a black hole with initial mass M(0), when accreting from 
a prograde accretion disk at a rate corresponding to 10% the Eddington ratio. Top pa nel: relat i ve ma ss 
increase, M{t)/M{0). Middle panel: spin parameter evolution, a, derived according to iBardeenI (Il970l) 's 
fomialism. Bottom panel: corresponding evolution of the radiative efficiency, e. 
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Fig. 3. — Face-on view of the inner part of the CND at t = 10 Myr. The color (logarithmic) scale refers 
to the gas surface density in units of 10^ Mq/pc^. Axes are in pc. Left panel: EVPle2-SFRTl-VHR-SS. 
Right panel: EVPle6-SFRT300-LowT. 
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Table 1 : Initial parameters for the different 
disk gas mass to lO^M©, and the stellar bulg 



runs: the black-hole mass is always fixed to 4 x 10 Mq, the 
e mass to 2 x lO^M©. 



Runs 


Particle number 


Particle mass [Mq] 


Evap. 


SF time 


Low-T 


RT 


Softening 




gas (bulge) 


gas 


(bulge) 


factor 


[Myr] 


cooling 




length [pc] 


EVPle2-SFRTl 


2 


X 10^ 


(10^) 


500 


(2 X 10'') 


lO'-^ 


1 


off 


off 


1.0 


EVPle2-SFRTl-LowT* 


2 


X 10^ 


(10^) 


500 


(2 X 10^) 


102 


1 


on 


off 


1.0 


EVPle2-SFRT300 


2 


X 10^ 


(10^) 


500 


(2 X lO'') 


102 


300 


off 


off 


1.0 


EVPle2-SFRT300-LowT* 


2 


X 10^ 


(10^) 


500 


(2 X 10*) 


102 


300 


on 


off 


1.0 


EVPle4-SFRTl 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


10" 


1 


off 


off 


1.0 


EVPle4-SFRTl-LowT* 


2 


X 10^ 


(105) 


500 


(2 X 10^) 


10" 


1 


on 


off 


1.0 


EVPle4-SFRT300 


2 


X 10^ 


(10^) 


500 


(2 X lO'') 


10" 


300 


off 


off 


1.0 


EVPle4-SFRT300-LowT* 


2 


X 10^ 


(10^) 


500 


(2 X 10*) 


10" 


300 


on 


off 


1.0 


EVPle6-SFRTl 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


106 


1 


off 


off 


1.0 


EVPle6-SFRTl-LowT* 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


10^ 


1 


on 


off 


1.0 


EVPle6-SFRT300 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


106 


300 


off 


off 


1.0 


EVPle6-SFRT300-LowT* 


2 


X 10^ 


(10^) 


500 


(2 X 10^) 


10^ 


300 


on 


off 


1.0 


EVPle2-SFRTl-IwO 


2 


X 10^ 


(10^) 


500 


(2 X 10'') 


102 


1 


off 


off 


1.0 


EVPle2-SFRT300-IW* 


2 


X 10^ 


(10^) 


500 


(2 X 10*) 


102 


300 


off 


off 


1.0 


EVPle2-SFRTl-LRt 


2 


X 10"^ 


(10^) 


5000 


(2 X 10^) 


102 


1 


off 


off 


1.0 


EVPle2-SFRTl-HR* 


2 


X 10^ 


(10^) 


50 


(2 X 10") 


102 


1 


off 


off 


1.0 


EVPle2-SFRTl-HR-SS^ 


2 


X 10^ 


(10^) 


50 


(2 X 10"*) 


102 


1 


off 


off 


0.1 


EVPle2-SFRTl-VHR-SS* 


2 


X 10^ 


(10^) 


5 


(2 X 10") 


102 


1 


off 


off 


0.1 


EVPle2-SFRTl-RT 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


102 


1 


off 


on 


1.0 


EVPle2-SFRT300-RT 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


102 


300 


off 


on 


1.0 


EVPle2-SFRT6 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


102 


6 


off 


off 


1.0 


EVPle2-SFRT10 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


102 


10 


off 


off 


1.0 


EVPle2-SFRT20 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


102 


20 


off 


off 


1.0 


EVPle2-SFRT50 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


102 


50 


off 


off 


1.0 


EVPle2-SFRT100 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


102 


100 


off 


off 


1.0 


EVPle2-SFRTl-SEt 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


102 


1 


on 


off 


1.0 


EVPle2-SFRT300-SEt 


2 


X 10^ 


(10^) 


500 


(2 X 10") 


102 


300 


on 


off 


1.0 



* The LowT label indicates runs with low-temperature metal cooling for Z = Zq. 

^ The IW label indicates runs with isoti'opic winds. The others are performed with axial winds. 
'^ The LR label stands for low resolution. 

* The HR label stands for high resolution; the HR-SS label stands for high resolution and small softening; 
the VHR-SS label stands for very high resolution and small softening (see text in Sect. l2.3l for more details). 
•f The SE label indicates runs with consistent stellar evolution calculations, and metal-dependent cooling in 
the temperature range ~ 10 — 10^ K, for Z evolving according to the proper stellar lifetimes and yields. 
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mass of 500 Mq. By comparison, the low-resolution (LR) run has the same gravitational softening 
length, but ten times less gas particles. The two high-resolution (HR) simulations have 2 million 
gas particles (50 Mq per particle) and are run with either a gravitational softening of 1 pc or with 
a smaller softening of 0.1 pc. Similarly, the very-high-resolution run (VHR) is run with a grav- 
itational softening of 0. 1 pc, but sampling our circum-nuclear disk with 20 million gas particles, 
corresponding to a gas mass resolution of 5 M©. 

To check more realistic values of star formation timescales, we re-run the previous standard cases 
also with t° = 6, 10, 20, 50, 100 Myr, besides the two t° = 1 Myr and t° = 300 Myr reference 
values. 

We then perform two simulations considering radiative transfer effects, as described in section IZTl 
with standard resolution, assuming Aq = 10^ and either t° = 1 Myr or t° = 300 Myr. 
Finally, we study the role of stellar evolution (SE) with consequent metal spreading by start- 
ing the simulations with a pristine gaseous environment (Z = 0), by polluting it with meta l 



yields (|Woosley & Weaver! Il995l: Ivan den Hoek & Groenewegen 



1997 



Thielemann et al.l l2003b. 



Padovani & Matteucci 



(119931) . and by 



from stellar evolution, according to the stellar lifetimes by 

consistently following the full chemical composition of the medium for H, He, O, C, Mg, S, Si, Fe. 

F or the technical details on thi s topic and studies on parameter dependencies, we refer the reader 



to 



Tomatore et al.l (12004 l2007l) : lMaio et al.l(l2007[l2010aLl201lh 



2.4. Star formation 



In all our runs, independently of the resolution, the softening, or other parameters, the disk experi- 
ences a burst of star formation, with a peak of a few solar masses per year, lasting for ~ 1 — 10 Myr 
(see figure H). These values are compatible, although on the low side, with the typical val- 
ues of star formation rates in circum-nuclear disks associated with star-bursts a nd late stages of 



galaxy interactions, from a few solar rnasses per year up to several hundreds (iKennicutt Il998 ; 



20041 : Evans et all 120061: lOownes & Eckart 



2007 ) 



Downes & Solomonlll998l : iDavies et al.ll2003 
After t ^ 100 Myr, the star formation rate decreases by 1-2 orders of magnitude, depending on 
the particular parameters used, but it always peaks at 3 — 5 Myr for t^ = 1 Myr, and the global 
behaviour is quite independent of Aq (left panel of figure ID and of resolution (right panel of fig- 
urelH). The cases with t° = 300 Myr result in a slightly later peak of star formation at ~ 5 — 30 Myr 
with slightly larger values. This is due to the longer timescale considered and the resulting higher 
dumpiness reached in the star forming regions. However, in all the simulations the central region 
of the disk is efficiently converted into stars, that become a relevant component after ~ 10 Myr. 
In figurelH we also show the star formation rates for different t° assumptions: t° = 1, 6, 10, 20, 50, 
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100, 300 Myr (left panel jj, and additional calculations including stellar evolution and metal pro- 
duction (right panel). We see (left panel of figure |5]) that when increasing t° - i.e. when lowering 
star formation efficiency - the SFRs during the initial phases (i.e., by < 10 Myr) get shallower and 
the bulk of formation of stars gets gradually shifted to later times. After the first burst, efficiently 
star forming models with short t° (say around < 50 Myr) show a clearly declining trend, due to 
the lack of remaining cold gas, that has been rapidly consumed during the initial phases. Feedback 
processes take place immediately, heat the surrounding medium, and inhibit further star formation. 
In less efficiently star forming models, with larger t° values, we assist to a slowly increasing trend 
of the SFRs for the initial tens Myrs. In these scenarios, gas can cool and condense for longer 
time without forming significant amounts of stars. At > 20 Myr, cooling processes finally lead to 
a number of strong SFR peaks, that take place when the medium is, in average, colder and more 
clumpy than in the previous cases. Thus, only after these episodes feedback effects become ef- 
ficient in heating the gas and halting further stellar production. Interestingly, a limiting situation 
is represented by t° = 50 Myr, for which the SFR is roughly flat over the entire lapse of time 
simulated, and reflects an almost continuous star formation process during which feedback effects 
and gas cooling are always well regulated. 

The inclusion of detailed stellar evolution with self-consistent metal production (right panel of 
figure |5]) determines some difference in the time evolution of the SFR due to the increasing metal- 
licities that slightly boost star formation. However, no major effects in the final results are found. 
In figure [6l we finally plot the evolution of the SFRs in models with different wind shapes (left) 
and different gas density thresholds (right). On the left panel, we show results from simulations 
with t° = 1 Myr and Aq = 10^, considering: axial winds with coupling with the ISM after a travel 
length of 1 pc (black solid line, simply labeled as "Axial winds"); axial winds with coupling with 
the ISM after a travel length of 10 pc (red dotted line, labeled as "Extended winds"); axial winds 
with coupling with the ISM after a travel length of 10 pc and density contrast below 0.1 (blue 
dashed line, labeled as "Extended low-density winds"); isotropic winds (cyan dot-dashed line). 
While the first and second case are basically very similar, some slight difference is found for the 
third case, where axial winds start interacting with the medium only in under-dense regions, and 
hence, lead to less star formation after the first few Myrs. In the case of isotropic winds, particles 
are not coUimated in one direction, but can interact with the surrounding gas at An, compressing 



'^ We note that a value of t^ = IMyr (with gas threshold of lO^'cm"'^) would correspond to star formation 
efficiencies of roughly r^ 50%, while <° — SOOMyr would correspond to a very low value of ~ 0.15%. As mentioned 
in the Introduction, these are quite extreme numbers, compared to expected values of a few per cents, or so. A more 
realistic assumption of star formation efficiency of ^ 1% would require i" ^ 50 Myr. Even though there are some 
cases (that we also take into account) for whi ch timescales can be significantly shorter, as in the Milky Way, for 



example, for which it is expecte d ^5 — 6 Myr (Burkertil2006), or the T-Tauri complex for which observed timescales 
are of the order of - 1 - 3 Myr jPalla & Stahlefcood iHartmann et al.lboOll) . 
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the gas in all directions, hence, a larger SFR is obtained. Quantitatively, the boost to the SFR due 
to isotropy is only a factor of a few, though. On the right panel, SFRs for the same initial setup and 
axial winds, but with gas density thresholds of 10^ — 10^ cm^^ are plotted. Since star formation is 
strongly dependent on density, larger thresholds determine higher SFRs and lower thresholds de- 
termine lower SFRs. In the cases considered, differences can reach up to ~ 2 orders of magnitude, 
ranging from < 10~^ Mq/ji (for a threshold of 10^ cm~^) up to several Mq/jic (for a threshold 
of 10^ cm^'^). The obvious consequence is that the little amount of stars formed in low-density- 
threshold scenarios will trigger feedback processes acting essentially on low-density gas, and that 
could be also effective in randomizing the dynamics of the medium (see more discussion in the 
parameter study later). 



3. Spin 

We take advantage of the extremely wide variety of our simulations to track the dynamics of the 
stream of particles inflowing towards the center and show the results for different cases. 



3.1. Method 

We estimate the accretion rate and angular momentum of the inflowing material through the small- 
est surface surrounding the SMBH that can be numerically resolveqj (e-g-, 3 pc and 0.3 pc for 
simulations with a 1 pc and 0.1 pc resolution, respectively). Due to the large range of scales in- 
volved, the formation of accretion disks around the SMBH could not be resolved down to ^ 0.01 
pc, where the gas can be modelled as a standard accretion disk. Therefore, the relevant proper- 
ties of the disk and their evolution in time have been simply inferred from the properties of the 
particles inflowing through the previously mentioned surface. In particular, we assume that: 1) 
the accretion rate onto the hole is equal to the accretion rate at the smallest resolved scales; 2) the 
direction of the orbital angular momentum in the external part of the unresolved accretion disk is 
parallel to the direction of the total angular momentum (relative to the SMBH) of the infalling gas 
at the resolution limit. We notice that all these assumptions are necessary working hypotheses. 
Nevertheless, we stress that these information can be gathered only by performing high-resolution 
numerical simulations that accurately resolve the "sphere of influence" of the black hole. We limit 
the gas inflow to the Eddington rate for the SMBHs, and use the SMBH accretion histories ob- 
tained from our SPH simulations to evolve the SMBH spin vector. We notice that while the SMBH 



We note that the sphere of influence (extending out to ^ 20 pc) is largely resolved in all our runs. 
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to 0.100 r 



LR 

HR 

HR 

VHR 



softening= 1 .Ope 
softenIng= 1 .Ope 
softening= 1 .Ope 
softening = 0.1 pc 
softening = 0.1 pc 




Fig. 4. — Star formation rates for simulations evaporation parameters, Aq (left), and resolution (right), as 
explained in the legends. The cases refer to t^ = 1 Myr with cooling down to ~ 10^ K. The initial drop is 
due to the denser central regions of the CND that undergo star formation and get evacuated during the first 
few Myrs. 



to 0.100 r 



5 

10 

20 

50 
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Myr 




tn O.'OO r 
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300 Myr, no SE 

1 Myr, SE 
300 Myr, SE 




Fig. 5. — Left panel: Star formation rates for simulations with different star formation time scales, t^ = 1> 
6, 10, 20, 50, 100, 300 Myr. Right panel: Star formation rates for simulations with inclusion of stellar 
evolution (SE) treatment and metal pollution process compared to the corresponding star formation rates in 
the runs with no SE treatment, for t° = 1 Myr and t° = 300 Myr, with ^o = 10^. 
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Fig. 6. — Left panel: Star formation rates for simulations with different wind parameters for the t^ = I Myr 
and Aq = 10^ case: axial winds with coupling with the ISM after a travel length of 1 pc (black solid line, 
simply labeled as "Axial winds"); axial winds with coupling with the ISM after a travel length of 10 pc (red 
dotted line, labeled as "Extended winds"); axial winds with coupling with the ISM after a travel length of 
10 pc and density contrast below 0.1 (blue dashed line, labeled as "Extended low-density winds"); isotropic 
winds (cyan dot-dashed line). Right panel: Star formation rates for simulations with gas density thi^esholds 
of 10^ — 10^ cm"'^, as indicated by the labels, for the t^ = 1 Myr and Aq = 10^ case with axial winds. 
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spin and mass evolve in the semianalytical model we use to postprocess our simulations (described 
in section [X2I ). the SMBH mass is constant in our SPH simulations, and the gas flowing toward 
the SMBH is not removed from the simulations. As shown in section [331 we limit our simulations 
to t < 20 Myr, to limit the mis-match between the SMBH mass in the simulations and in the 
semianalytic model to a factor of a few. We stress that such a small mass growth would not change 
substantially the gas dynamics, since the CND is significantly more massive than the SMBH in our 
initial set-up. 
If the orbital angular momentum of the disk around the SMBH is misaligned with respect to 



the SMBH spin, the coupled acti on of viscosity and re 



the disk in its innermost re gion (IBardeen & 



viscous/accretion timescale (Scheuer & Feiler 



'etterson 



ativist ic Lense-Thirring precession warps 



1975) over a timescale shorter than the 



1996t Natarajan & Pringlelll998l : lBogdanovic et al 



20071). We develop an analytical framework to model the response of the accretion disk (in the 
innermost non-resolved region) to the angular momentum of the mass inflow as described below. 
At each timestep v^e can calculate th e evolution of the magnitude and direction of the SMBH spin 
JPerego et al.ll2009l : botti et alUlQlol) . 



3.2. Bardeen-Petterson effect 



We use the SMBH accretion histories obtained from our SPH simulations to follow the evolution 
of each SMBH spin vector, Jbh = (clG M^yl/'^) Jbh, where < a < 1 is the dimensionless spin 
parameter and jbh is the sp in unit vector. We will use the term maximally rotating for a black hole 
with a = 0.998, following iThomel (|l974l) who showed that accretion driven spin-up is limited to 
such value. Magnetic fields connecting material in the disk and the plunging region may further 
reduce the equilibrium spin. Magneto-hydrodynamic si mulations for a serie s of thick accretion 
disks suggest an asymptotic equilibrium spin at a ~ 0.9 (IGammie et al.ll2004l) . 



The sc heme we adopt to study spin evolution is based on the model developed in IPerego et al 



(120091) . Here we summari ze the algorithm. We assume tha t inflowing gas forms a geometrically 
thin/optically thick «-disk (|Shakura & Sunyaevlll973l Il976|) on sub-parsec scales (not resolved in 



the simulation), and that the outer disk orientation is defined by the unit vector Icdgc- The evolution 
of the a-disk is related to the radial viscosity, z/i, and the vertical viscosity, 1^2'- ^i is the radial 
shear viscosity, while U2 is the vertical shear viscosity associated to the diffusion of vertical warps 
through the disk. The two viscosities can be described in terms of two different dimensionless 
viscosity parameters, ai and 02, through the relations 



i^l,2 = ttl,2^Cs 



(7) 



where H is the disk vertical scale height and Cg is the sound speed of t he gas in the accretion disk . 
Both of them depend on R, Mbh, M and ai. Here we assume the IShakura & SunyaevI (|1973h 
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solutions for the ga s pressure dominated zon e. 
We further assume ([Lodato & Pringlell2007|) 



with@ 

and /2 = 0.6, i.e. 



«2 = /2/(2ai), 

ai = 0.1 
«2 = 3, 



and power law profiles for the two viscosities (jShakura & Sunyaevlll973Lll976r) . 

Z/1.2 OC i?3/4. 



(8) 

(9) 
(10) 

(11) 



If the orbital angular momentum of the disk around the SMBH is misaligned with respect to the 
SMBH spin, the coupled action of viscosity and relativistic Lense-Thirring precession warps the 
disk in its innermost region f orcing the fluid to rotate in the equatorial plane of the spinning SMBH 
(IBardeen & Pettersonlll975l: Kumar & Pringlelll985|) . The timescale of propagation of the warp is 



short compared with the viscous/accretion timescale (IScheuer & Feilerll 19961 Natarajan & Pringle 



19981), so that the deformed disk reaches an equilibrium profile that can be computed by solving 
the equation of conservation of angular momentum 



RdR^ ^ RdR \ dR 



1 d fl dl 

RdR [r'^^ OR 



+ 
2GJbh X L 



i?3 



(12) 



where vji is the radial drift velocity, S is the surface density, fi is the Keplerian velocity of the 

gas in the disk, and L is the local angular momentum surface density of the disk, defined by its 

modulus L and the versor 1 that defines its direction. 

The boundary conditions to eq. (fT2]) are the direction of L at the outer edge ledge, the mass accretion 

rate (that fixes the magnitude of S), and the values of mass and spin of the massive BH. 

The values of ledge and of the mass accretion rate are directly extracted from the simulations: 

in particular, they are computed considering those SPH particles nearing the MBH gravitational 

sphere of influence that are accreted. 

The spin value, Jbh, is treated on post-processing in a Monte Carlo fashion, by randomly assuming 

several different realizations at the initial time, and evolving each of them. More exactly, the 

direction of the SMBH spin changes in response to its gravito-magnetic interaction with the disk 



' These values are slightly different from others also found in hterature (as e.g. in lBate et al.ll2010l) 
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on a timescale longer than the time scale of warp propagation (|Perego et al.ll2009|) . This interaction 
tends to reduce the degree of misalignr nent between the disk an d the SMBH spin, decreasing with 
time the angle between Jbh and ledge (|Bogdanovic et alj|2007|) . So, the SMBH spin evolution is 
followed by solving for the equation 

dJBH ,Va/7-. M/7-. \ ^ttG f L X Jbh 



dt 



MA(i?iso)l(i?iso) 



disk 



/22 



-dR. 



(13) 



The first term in eq. (fT3T ) accounts for the angular momentum deposited onto the SMBH by the 
accreted particles at the innermost stable orbit (ISO), where A(i?iso) denotes the specific angular 
momentum at -Riso and l(i?iso) the unit vector parallel to Jbh^ describing the warped shape of 
the disk. The second term instead accounts for the gravo-magnetic interaction of the SMBH spin 
with the warped disk. It modifies only the SMBH spin direction (and not its modulus)^ conserving 
the total angular momentum of the composite (MBH+disk) system (|King et al.ll2005h . The inte- 
grand peaks at the warp radius (-Rwarp). where the disk deformation is the largesllj. The equation 
incorporates two timescales: the accretion time related to the first right-hand term describing the 
e— folding increase of the spin modulus, and the shorter timescale of SMBH spin alignment 

-2/35 

/Edf'Vr, (14) 



T-ai ~ lO^a' 



5.5/7 



M, 



BH 



^4 X lO^M©^ 

that will ensure a high degree of SMBH-disk gravito-magnetic coupling during SMBH inspiral. In 
Eq (fT4l) . /Edd is the SMBH luminosity in units of L^dd- 

We apply iteratively eq. (fT2l) and (fT3T ). by using inputs from the SPH simulations that give the 
values of the mass accretion rate, the SMBH mass and the direction of ledge • The algorithm returns, 
as output, the spin vector, that is, its magnitude and direction. At each timestep our code therefore 
provides the angle between the spin vector of the SMBH and the global angular momentum vector 
of the circum-nuclear disk. 



3.3. Resulting spin alignment and magnitude 

We will now discuss the results expected for the spin behaviour in the various runs presented be- 
fore. We will start with general trends and then show also one of the highest resolution runs. 
Figure |7] shows the time evolution of the spin magnitudes (top panel). We show here the results 
of 100 Monte Carlo realizations starting from a flat distribution in both spin magnitudes and ini- 
tial orientations for four representative cases (the alignment process and the evolution of the spin 



^ The exact definition of i?warp is where the vertical viscous timescale t^.-^ ~ R^ 1^2 in the disk is comparable to 
the Lense-T hirring precession t imescale. Because i?warp and the radius at which the disk is maximally deformed are 
comparable (IPerego et al .112009b . we simplify the notation in the paper using only R^ 



^warp ■ 
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Fig. 7. — For each panel, the time evolution of the spin parameter (top), of the spin orientation angle 
(middle), and of the luminosity in units of the Eddington luminosity, f^dd (bottom), computed at a distance 
of 3 pc from the central BH, are shown. In the upper and middle panels the solid lines refer to the mean 
values of the spin magnitude and orientation 6, averaged over 100 realizations. When 9 approaches zero, 
the accreted material is aligned with the SMBH spin ("coherent" accretion). When 6 is lai^ge, material is 
flowing with angular momentum misaligned with respect to the spin, causing the spin to decrease. The 
darker and lighter (red and orange in the online version) shaded areas refer to the regions encompassing 
1- and 2-sigma around the means. The different panels refer to the runs with 500 Mq mass resolution: 
EVPle2-SFRTl (top left), EVPle2-SFRTl-LowT (top right), EVPle2-SFRTl-IW (bottom left), EVPle2- 
SFRTl-RT (bottom right). In the top right panel the evolution of f^dd for the EVPle2-SFRTl run is shown 
with a dashed Une to facilitate the comparison. 
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parameter are almost identical in all other runs). The time evolution of the relative angle, 9, be- 
tween the spin of the SMBH and the orbital angular momentum of the circum-nuclear disk is also 
shown (central panel). The inner (outer) shaded areas in the upper and central insets correspond 
to 1-sigma (2-sigma) deviations from the mean values, over 100 random initial spin realizationsl^. 
The corresponding luminosity in units of the Eddington luminosity, fEdd, is plotted, as well (bot- 
tom panel). All the results discussed in the following, where not specified otherwise, are obtained 
computing the properties of the accreting material at a distance from the central SMBH of 3 pc, i.e. 
three softening lengths of the low resolution runs, and for star formation timescale of t° = 1 Myr. 
Top-left panels refer to the reference run with axial wind and radiative cooling down to ~ iC K, 
top-right panels to the corresponding run with low-temperature cooling (below 10'' K), bottom-left 
panels to the run with isotropic winds and cooling down to ~ 10'' K, and bottom-right panels to 
the reference run including RT. 

In general, SMBH spins rapidly loose memory of their initial orientation, and accretion torques suf- 
fice to align the spins with the angular momentum of the disk orbit on a short timescale (;$ 1 — 2 
Myr). In all the cases, Jbh tends to align with the angular momentum of the accreting gas, as 
extracted by the simulations and highlighted by the black lines in the upper and central insets, 
within a few Myr, regardless the initial orientation and magnitude of the SMBH spin, even if ini- 
tially completely anti-aligned. After the first, fast alignment of the spin with the orbital angular 
momentum of the circum-nuclear disk, the accretion flow is mostly coherent, i.e. the inflow traces 
the angular momentum of the large-scale circum-nuclear disk that feeds the SMBH, and the gas 
luminosity near the center is very close to the Eddington luminosity. Only in the RT runs (bottom- 
right panel) we find inclination angles that are slightly larger: this is due to the more realistic 
treatment of entropy injection by radiative sources, that is powered by radiation heating (missing 
in the other cases). This does not change our overall conclusions, because typical values are still 
< 50°, and only rarely isolated spikes are found. After the "fast" spin alignment, accretion pro- 
ceeds essentially coherently, i.e. the gas accretes on planes that have an inclination angle with 
respect to Jbh that is significantly smaller than 90°. As a consequence, a monotonically increases 
with time, except for a very close to unity, where it stays roughly con stant, or decreases slightly. 



Since for a ~ the alignment timescale (oc a^/^) is extremely short (|Perego et al.ll2009h . slowly 
rotating SMBHs start immediately to accrete in a prograde fashion, and the lower boundary of the 
shaded areas are not expected to decrease with time on timescales of ~ Myr. Once a BH starts 
accreting from a prograde disk, the spin magnitude also increases over short timescales, close to 
the maximum allowed one. 

The resulting accretion rate, fEdd, shown in the lower insets, is more seriously affected by the 
different physical mechanisms By comparing the EVPle2-SFRTl and the EVPle2-SFRTl-LowT 



We checked that changing the number of initial realizations the final results are not affected significantly. 
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runs (upper right panel) it emerges that low-temperature cooling contributes to a slight enhance- 
ment of /Edd- More quantitatively, the average increment during the first 15 Myr is roughly ~ 4 
per cent, with temporary peaks of a factor of ~ 2. In fact, at ~ 5 Myr, f^dd results to be ~ 0.3 
for the standard run, and ~ 0.6 for the low-temperature case. This is related to the more effi- 
cient gas condensation process that can take place in such conditions. In this case, gas can loose 
more pressure support, reach lower temperatures, and form smaller clumps with very high densi- 
ties on short timescales. Both formation of dense structures and gas compression from SN shocks 
reshuffle the local angular momentum of the disk strengthening inflows onto the central black hole 
Jwadaetal.ll2009l : iKawakatu & Wadall2009l : iHobbs et allboi lb . The overall effect is limited on 
short timescales, since in the long run (e.g. > 15 Myr) the trends tend to converge mainly because 
of gas consumption, as shown by the decrement in the star formation activity. 
When comparing isotropic winds to axial winds (see upper-left and lower-left panels), we notice an 
initial drop in accretion rate for the isotropic case. The very first episodes of star formation happen 
close to the black hole, where the gas density reaches initially its peak. Because of the disk-like 
geometry of our gas distribution, isotropic winds can interact with close gas particle more effi- 
ciently, initially evacuating the central region of the disk, and hence decreasing the accretion rate. 
However, as soon as star formation picks up in the outer regions of the disk, the accretion rate onto 
the black hole raises again, up to values even slightly higher than in the runs with axial winds. 
This is because isotropic winds mix the angular momentum of gas more efficiently. In the RT case 
(bottom-right panel), heating from radiative feedback keeps the gas hotter, strongly suppressing 
the formation of dense substructures and star formation. As a consequence, /sdd is decreased by a 
factor of a few, with peak values of f^dd ^1/2- 

In our highest resolution run (figure |8]), we can resolve almost one star at a time. Even in this run 
we find results similar to the lower resolution runs, with the only difference that we can follow gas 
inflow down to sub-parsec scales (0.3 pc), deeper in the potential well of the central SMBH, where 
the accretion rate is, in average, a factor of ~ 2 — 3 larger. 

We emphasize that these results apply only to quasars in merger remnants, not to those fed by 
disk instabilities and cold streams. In tho se cases, massive clumps may chaoticize accretion (e.g. 
Boumaud et al.ll201ll : Dubois et al.ll2012l) . since they would usually infall from large scales carry- 
ing their own angular momentum. 



3.4. Parameter study 



To draw more robust and definitive conclusions, we finish by showing the main findings arising 
from a parameter study for different cases. In particular, to check the dependencies on star forma- 
tion we explore the changes induced by adopting different star formation timescales, gas density 
thresholds, evaporation factors, self-consistent stellar evolution, and metal spreading. 
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Fig. 8. — Same as figure |7]for the highest resolution run (simulation EVPle2-SFRTl-VHR-SS). 




Fig. 9. — As in figure|7]for runs performed assuming different star formation timescales: t^ = 50 (EVPle2- 
SFRT50, left panel), 100 (EVPle2-SFRT100, middle panel), and 300 Myr (EVPle2-SFRT300, right panel). 
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Similarly to figure |7] and |8l where we explored the case t° = 1 Myr, in figure |9] we display the 
time evolution of the spin parameter, of the spin orientation angle, and of the luminosity in Ed- 
dington units, /Edd> also for the runs performed with t° = 50, 100, 300 Myr (in the cases with 
t° = 1 — 20 Myr no relevant differences are found, hence we will not focus on each of them). The 
general trends show that, despite the good agreements of the spin parameters (in all the cases rang- 
ing between ~ 0.6 and 0.9 at t ~ 10 Myr) and of the inclination angles (always quite small), mod- 
els with larger star formation timescales (i.e. smaller efficiency), as the one with t° = 300 Myr, can 
nevertheless accrete slightly more mass and sustain f^dd at slightly larger values (/Edd ~ 0.8 — 1 
at > 15 Myr) than models with more efficient star formation and shorter t°. We also notice that, 
in models with larger t°, the formation of long-lived gas clumps that could stir the medium seems 
not to be very efficient in increasing the chaotic behaviour of the CND, with /sdd reaching roughly 
unity after > 10 - 15 Myr. 

Our conclusions could be affected by the ability of resolving properly the cold, dense, star forming 
regions. As a check, we analyze the cases carried out by assuming different gas density thresholds 
for star formation in the range 10^ — 10^ cm~^ (see supplementary online material). It emerges 
that models with low-density thresholds have a more chaotic behaviour and show larger spreads in 
the results. Indeed, for the scenarios with a threshold of ~ 10^ — 10'^ cm"'^, the spin parameter is 
roughly a ~ 0.4 — 0.95 and 9 can smoothly increase above ~ 30°, or so. The more chaotic behav- 
ior is visible mainly at initial times, when first bursts of star formation take place, and the spreads 
in the data are wider than in the higher-density-threshold cases. Thus, the accretion rate is eas- 
ily inhibited (with /gdd dropping below ~ 0.1), because even material at relatively low-densities 
experiences star formation and gets evacuated by the consequent thermal heating and feedback 
processes. Larger, more realistic, density thresholds, instead, allow the gas of the CND to cool 
and actually condense, so that star formation takes place only in very clumpy regions on the disk. 
Dense, cold fragments are affected in a minor way by feedback processes, can survive them, and 
continue accreting on the central BH. As a consequence, the CND can keep a more stable and 
ordered shape, and the cases with a threshold > 10^ cm""^ converge rapidly (in < 1 Myr) to similar 
trends. Moreover, they are compatible with the fiducial value of ~ lO'' cm~^ used throughout this 
paper. 

The effects of changes in the evaporation factor (see supplementary online material) lead to sim- 
ilar spin evolution, irrespectively from the detailed values used, and to inclination angles always 
below ~ 50°. In general, despite the non- very-regular trends, some boosts in 9 are present when 
there is strong star formation activity (figure 5) for all the cases, with corresponding decrements in 
the Eddington fraction, since hot gas gets more pressure supported and provokes smaller inflows 
towards the central regions. Instead, when star formation rates drop down, the resulting 9 values 
decrease and the /sdd values get closer to unity. This variability is a consequence of the interplay 
between star formation heating and gas cooling that regulate the behaviour of the environment. 
Stellar feedback heats the medium, while simultaneous cooling processes cool down the gas, sta- 
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bilize the disk, and trigger the next star formation episodes. 

Complete full stellar evolution and metal spreading, displayed in figure[10l does not imply dramatic 
changes in the spin results. The sharper variations (on timescales of ~ 1 Myr) of the inclination an- 
gle and of /Edd are due to the larger radiative rates in presence of metals, and, hence, faster cooling 
and star formation processes, that respectively boost and inhibit gas inflows. This is mostly evident 
for t1 = 1 Myr, where metal pollution takes place earlier than in the case of longer t° = 300 Myr. 
In the former case, the effects of short star formation timescales (see also right panel in figure |5) 
are enhanced by the stronger radiative losses of the newly enriched material, that quickly cool the 
ambient medium and temporarily enhance gas infall. In the latter case, the star formation rate at the 
same times is a factor of a few to ~ 10 more limited, and cooling from metals is able to efficiently 
stabilize the disk against the chaotic motions from the rarer stellar feedback. These conclusions on 
the stabilizing role of metals are easily cross-checked by a direct comparison with the correspond- 
ing Z = cases, in the top-left plot of figure |7] and in the right plot of figure |9l for t° = 1 Myr 
and t° = 300 Myr, respectively. 



4. Discussion and conclusions 

We have presented a number of thre-dimensional, N-body, hydrodynamical, chemistry simulations 
of the behaviour of circum-nuclear discs around super-massive black holes. We considered the 
relevant local physical processes that could alter the coherency of gaseous flows onto the central 
regions, i.e.: cooling, star formation, feedback effects, and radiative transfer from stellar sources. 
Our goal was to explore the implications of these several phenomena on the stability of circum- 
nuclear discs around SMBHs and the consequent effects on the black-hole spin evolution. Our 
investigation considers mainly those "local" processes that could increase the turbulence of the 
gas (such as the gas self-gravity, cooling, star-formation, and feedback from young stars), while 
global processes, such as the precession of the nuclear gas structures in a less symmetric, tri-axial 
potential, are not taken in account, as we started with an axisymmetric distribution of gas and stars. 
Global instabilities or asymmetries that c ould decrease the degree of coherency in the BH fuelling 



(as discussed e.g. in lHopkins et al.ll2011l) are beyond the aims of this work. 
Our study does not include AGN feedback. The energy liberated by an AGN may be by far larger 
than the typical energies delivered through stellar feedback in the central regions of the galaxy. 
The effectiveness of the AGN in altering the CND gas dynamics, however, depends on the jet ge- 
ometry and on the fraction of energy actually deposited in the nuclear regions. In general terms, 
we can make some inferences on the expected impa ct of AGN feedback from the results of (two- 



dimensional) dedicated simulations such as those by iNovak et al.l (|201 1|) . They find that mechani- 
cal feedback (winds) dominates the effect on the local gas dynamics, and that this effect is stronger 
when the wind is confined in a cone. Nevertheless, even during a high-accretion event a beamed 
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wind does not strongly affect the gas located perpendicularly to the outflow. These simulations 
differ from our three-dimensional setup, though. On the one hand, they do not include any dense 
and rotationally supported gas structure, that would reduce even further the effect of the AGN 
feedback. On the other hand, they keep the orientation of the conical outflow constant. If the 
orientation would follow the rearrangement of the SMBH spin, it would in some cases transiently 
impinge on the CND, modifying more strongly its dynamics and the properties of the accreting 
flow. Definitely, more self-consistent investigations are needed to assess the interplay between 
AGN feedback and CND dynamics, but this first study, focusing on the dynamical perturbations 
self-generated by the CND, will facilitate distinguishing the effects of AGN feedback in future 
works. 

Summarizing, our analysis suggests that in most of the cases local star formation episodes and 
feedback effects are able to partially break the general coherency of the gas flow, but do not inject 
enough energy and momentum into the gas near SMBHs to create a strongly chaotic environment 
that completely randomizes the orbits of the inflowing particles. As a conclusion, we find that 
initially maximally rotating SMBHs are slightly spun down, and initially slow-rotating black holes 
are spun up, leading to upper-intermediate equilibrium values, a ~ 0.6 — 0.9 (corresponding to 
radiative efficiencies e ^ 9% — 15%). 

Considering the whole parameter space, it results that different feedback mechanisms can affect 
the CND gas dynamics and coherency in different ways. Mechanical (kinetic and thermal) feed- 
back processes are responsible for causing temporary, moderate, changes in the inclination angle 
and Eddington fraction. Metal feedback, instead, introduces sharper variations, but has a generally 
stabilizing role for the CND gas flows, since it facilitates the formation of cold clumps. Radiative 
feedback is responsible for inducing a more chaotic regime with increments of the alignment angle 
of at least a factor of a few (up to ~ 50°), and for limiting the accretion capabilities of the central 
BH of a factor > 2. 

We note that our results are based on the angular momentum of the material that is effectively 
accreted on the SMBH, and not on the global ali gnment of the nuclear disc with the large scale 



properties of the host galaxy (|Hopkins et al.ll201 1|) . While the nuclear disc can be misaligned with 



respect to the galaxy disc, the relevant quantity to address spin evolution is the direction of the 
angular momentum of the material within the nuclear disc that eventually feeds the SMBH, not 
the alignment or misalignment between nuclear disc and the host. These results are in line with 



the su ggestion that jets are preferentially perpendicular to dust lanes (|Verdoes Kleijn & de Zeeuw 



20051), which are believed to be the remnants of galaxy mergers - a low-mass analogue of nuclear 



discs. Further hints for at least partial alignme nt are obtained fro m analyses of active galactic nu 



clei (AGN) and their hosts in the Sloan survey ([Lagos et al.ll201ll) . According to these studies, the 



material feeding the black hole is expected to accrete in a somewhat coherent fashion. 
Observationally, spin magnitudes of SMBHs powering AGNs can be estimated directly by mea- 
suring the width of the Kq, iron line, at 6.4 keV, through X-ray spectroscopy: recent determinations 
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with Chandrcii XMM-Newtoti^, and Suzaki\!_!\ ha ve revealed variable, re lativistic i 



iron emission 



lines from the inner disc in many Seyfert galaxies ([Brenneman et alj|201 l|i 'H It is suggestive that 



all SMBHs w ith spin measured directl y through X-ray spectroscopy have spin parameters between 



~ 0.6 and 1 (|Brenneman et alj|201 Ih . in broad agreement with our results. Indeed, if after the 
quasar phase the mass of a SMBH has not changed much through gas accretion or SMBH mergers, 
the expectation is that today's spins are similar to the final spin gained during the last growth spurt. 
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Fig. 10. — As in figure [7] for runs assuming star formation timescales of 1 (left) and 300 (right) Myr, 
both with full stellar evolution and metal enrichment calculations (SE). For a direct comparison with the 
corresponding Z = cases, see top-left plot in figure |2l and right plot in figure |9] respectively. 



